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Abstract: We demonstrate how to efficiently implement extremely 

high-dimensional compressive imaging of a bi-photon probability dis¬ 
tribution. Our method uses fast-Hadamard-transform Kronecker-based 
compressive sensing to acquire the joint space distribution. We list, in de¬ 
tail, the operations necessary to enable fast-transform-based matrix-vector 
operations in the Joint space to reconstruct a 16.8 million-dimensional 
image in less than 10 minutes. Within a subspace of that image exists a 
3.2 million-dimensional bi-photon probability distribution. In addition, we 
demonstrate how the marginal distributions can aid in the accuracy of joint 
space distribution reconstructions. 
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1. Introduction 


Characterizing high-dimensional joint systems is a difficult problem due to experimental im- 
practicalities such as long measurement times, low flux, or insufficient computing resources. 
One example of such a characterization is of continuous-variable entangled states - a resource 
gaining ground in quantum technologies . A widely used source of continuous-variable en¬ 
tangled states is Spontaneous Parametric Down-Conversion (SPDC) in a nonlinear crystal Q. 
Depending on the configuration, the resulting bi-photon state may be entangled in the transverse 
degrees of freedom To determine if the system is entangled, both the bi-photon joint 

position and joint momentum probability distributions composed of signal and idler photons 
must be measured through correlation measurements. 

Much work has been done recently to characterize high-dimensional position-momentum en¬ 
tanglement with discrete measurements l IspTO . Characterizations of position or momentum 
distributions resulting from SPDC are done by measuring signal and idler pixel correlations 
in either an image plane of the crystal (constituting a position measurement) or a Fourier- 
transform plane of the crystal (constituting a momentum measurement) through coincidence 
counting. 

For a bi-photon probability distribution, measurements are typically done by measuring cor¬ 
relations via raster scanning through the individual signal- and idler-photon probability dis¬ 
tributions. The time required to complete a raster scan with single-photon detectors quickly be¬ 
comes impractical for certain scans. Imaging these distributions with a camera has been shown 
in |T3[Tg, yet cameras often introduce far more noise than a single-photon detector. 

Recently, Compressive Sensing (CS) techniques have been introduced as an alter¬ 

native to raster scanning for characterizing a high-dimensional entangled system, dropping the 
measurement time from months to hours | 2T|[22) . While the data-acquisition time is drastically 
reduced, it comes at the cost of computational complexity, requiring a computational recon¬ 
struction of the signal. Performing CS on high-dimensional signals is not a new problem, and 
several clever solutions exist for utilizing separable compressive sensing matrices combined by 
a Kronecker product |23 However, these methods are ill suited for sampling the correla¬ 
tions in a joint space. 

In this article we propose the use of fast Hadamard transforms for high-dimensional joint 
space reconstructions. Specifically, we show how the Kronecker-product-based recursion rela¬ 
tions of Sylvester-type Hadamard matrices can combine single-particle sensing matrices. This, 
in turn, enables the use of fast Hadamard transforms in the joint space as they have been shown 
to drastically reduce CS reconstruction times | |25) . Using the randomization techniques outlined 
in | |26|[27) , sensing matrices composed of randomized Hadamard matrices offer tremendous 
speed enhancements in many reconstruction algorithms. 

Additionally we show that by using the individual signal and idler marginal distributions in 
our reconstruction of the joint space distribution, we can more accurately acquire transverse 
spatial correlations as measured by the mutual information. To the best of our knowledge, this 
is the first time the single particle information has been used in reconstructing the joint space 
distribution. 

To demonstrate the effectiveness of our method, reconstruct a 16.8 x 10® dimensional joint- 
space distribution in only a few minutes. Within this joint-space lives a 3.2 million-dimensional 
bi-photon position probability distribution from which we measure the degree of transverse 
correlations between signal and idler photons using the mutual information. 

The experimental realization here is closely related to the work performed in pT) . The exper¬ 
iment in this article is merely meant to demonstrate how structured randomness enables efficient 
reconstructions of the joint space distribution at even higher dimensions. In theory, this increase 
in resolution allows for an increase in the amount of measurable mutual information. 













2. Compressive sensing 

CS helps to overcome unreasonable data-acquisition times associated with sampling signals 
with limited resources. CS requires that the signal of interest x has a sparse representation x 
via a basis transform. Limiting the discussion to real signals for simplicity, x G then there 
exists a basis transform 'L in which k <N components of x = 'L • x are nonzero; x is a A:-sparse 
representation of x. CS posits that if x is approximately A:-sparse, then only M =0 {k\og{N/k)) 
projections of x are needed to accurately sample x. A typical CS technique involves taking M ^ 
N projections of a signal x with a random sensing matrix A G to form a measurement 

vector y = A • X where y G Assuming the signal representation x is k-sparse, x may be 
recovered by solving 

min Tg(x) subjectto lly —A--x) ll? < e, (1) 

where the second term is a least-squares term bounded by a predefined error £ and 

Np:= 

is defined as the Zp-norm. The function g{x) is a function to be minimized depending on the 
assumed sparsity; g{x) is often chosen to be the /i-norm ||x|| j, the total variation as defined by 
the gradient || Vx|| i, or a combination of functions where x is known to be sparse. T is simply 
a weighting term that weights the sparsity, favoring either the least-squares solution or the 
sparsity. Because A G where M <^N, there are an infinite number of solutions confined 

to the least-squares term. The function g{x) picks out the sparsest of these solutions which 
should correspond to our signal x. An overview of compressive sensing and its applications 
may be found in p8) . 


3. Measuring a non-separable joint system 

We apply CS to measure the joint position probability distribution of the down-converted signal 
and idler photons from SPDC. Quantum mechanics tells us that the bi-photon state exists in a 
Hilbert space composed of the tensor product of the individual signal and idler photon Hilbert 
spaces. After representing the bi-photon state in a discrete basis, operators on these states are 
represented as matrices. In order make a measurement, we approximate the state as living in 
a finite dimensional Hilbert space. We can therefore represent a bi-photon operator matrix in 
terms of Kronecker products of individual signal and idler photon operator matrices. For CS, 
we can let these operators be projection operators and manipulate them such that they form the 
rows of a sensing matrix A. The sets of projection operators for the signal and idler spaces are 
designated by a a set of M patterns where each pattern is in Pj and P/ G In this 

manner, the sensing matrix is written as 


P5[1]®P/[1] 

P5[2]®P/[2] 


Ps[M](g)P/[M] 


( 2 ) 


for i G 1...M where P[/] represents the row of P. 

The Kronecker product ® in this article operates on matrices and vectors such that if a is of 
dimension mx n and b is of dimension p x q, then their Kronecker is of dimension mp x nq 





Fig. 1. The above experimental diagram demonstrates how to image a joint two-particle 
system. In this paper, the joint system is composed of highly-correlated signal and idler 
photons from a SPDC source. The experiment samples the position distribution of the joint 
system by taking random projections of signal and idler intensities with a spatial light 
modulator within an image plane of the crystal. An avalanche photodiode (APD) detects 
photon arrivals while the photon counters measure photon coincidences. 


represented as 



Q.\n\i 

Clm\^ 

Qmn}^ 


An experimental diagram for measuring the joint space bi-photon position probability dis¬ 
tribution is presented in Fig.[2where each particle’s space is depicted as two-dimensional. Yet, 
to simplify the CS formalism, we represent the signal and idler spaces as one-dimensional liv¬ 
ing in and the joint space distribution as a vector x G . As outlined in JiT) , compressive 
sensing is experimentally accomplished by taking random projections with patterns composed 
of [■\/N X '/n] pixels within each subspace of the signal and idler systems and then measuring 
the resulting correlations in photon counts. 

During reconstruction, the bi-photon probability distribution is already sparse in the pixel 
basis due to the tight pixel correlations resulting from energy and momentum conservation. 
This eliminates the need to define a sparse-basis such that 'T becomes the identity operator 
‘T = 1 in the formalism above. 

Random binary matrices were used in pTj , yet we wish to use structured random binary ma¬ 
trices for reconstruction purposes. Using properties of Kronecker products enables relatively 
efficient computations of the reconstruction operations A • x and • y, where A^ is the trans¬ 
pose of A, because A never needs to be computed explicitly p9| . However, structured random¬ 
ness can enable the use of fast transforms which are even more efficient. Our contribution is to 






demonstrate in the following section how randomized Sylvester-Hadamard matrices enable the 
use of fast Hadamard transforms in joint space CS reconstructions. 

4. Fast Hadamard transform based sensing matrices 

4.1. Randomly sampled & permuted Hadamard sensing matrices 

Sylvester-Hadamard matrices have a structure that is particularly advantageous to the CS frame¬ 
work. These matrices are generated from a simple recursion relation defined by a Kronecker 
product. 


Hi 

H 2 


[ 1 ] 

■ 1 1 

1 -1 


From these, any Sylvester-Hadamard matrix can be decomposed as follows; 


H 2 A — H 2 H2k— 1 


H 2 /:-1 H2A:— 1 

H 2/;-1 


(4) 


for k > 1. Because of this structure, Sylvester-Hadamard matrices are restricted to powers of 
two but can be used to build patterns and a sensing matrix that utilizes the speed and efficiency 
of a fast Hadamard transform 1K[...]. We use a normally-ordered fast transform in this work. Its 
algorithm is similar to that of a fast Fourier transform, but it consists of only additions and sub¬ 
tractions. Hence, it performs the reconstruction operations A • x and • y in 0 (N^ logA) time 
- significantly faster than an explicit matrix-vector multiplication of 0 (N^) time. A thorough 
overview of Hadamard matrices, fast-Hadamard transforms, and their applications to signal and 
image processing may be found in 

To construct P 5 , P/, and A from Hadamard matrices, the Hadamard matrices must be ran¬ 
domized in both their rows and columns. The sensing matrices must be both incoherent with 
the image yet span the space in which the signal resides. In other words, the sensing matrices 
must adequately sample the basis components of the signal. Random sensing matrices perform 
this task well, yet Hadamard matrices naturally contain much structure. To begin, each sensing 
matrix must be formed by taking specific rows from a Hadamard matrix with the correct di¬ 
mensions. A is constructed from H ^2 while P 5 and P/ are constructed from H^r. Because of the 
relation in Eq. (|^, the rows of A will be determined by the rows of P 5 and P/. 

The randomization of the Hadamard matrix rows is accomplished by defining two vectors 
rj and r/ G for each signal and idler system composed of M randomly chosen integers 
on the interval [2,N]. The values in r state which rows should be extracted from H^r when 
constructing P 5 and P/. Note that the interval begins at 2 because the first row of a Hadamard 
matrix is composed entirely of ones. The interval may begin at 1 if the total photon flux on 
a detector is desired. Also, note that A G where M << N^. This condition allows for 

scenarios where Ps,P/ G K'ttxiV M > N in the individual subspaces, meaning rows of 

Ha^ may be repeated when constructing P 5 and P/. 

The randomization of the Hadamard columns is accomplished by defining permutation vec¬ 
tors p 5 and p/ G that randomly permute the N columns of H^r. Once r and p have been 
defined for both the signal and idler subspaces, patterns are constructed by the following equa¬ 
tions: 


Ps = 
P/ = 


Hiv[rs,p.y] 

HAr[r/,P/] 


( 5 ) 






where the y and x components of H[y,x] refer to the rows and columns of H respectively. 

Although these operations are defined for the signal and idler subspaces, they combine in 
a particular way according to Eq. Q to construct a Hadamard-based sensing matrix A that 
enables fast transform operations in the joint space. The manner in which they combine to 
manipulate a Hadamard matrix H ^2 that spans the joint space is detailed in the next section. 

4.2. Joint space Sylvester-Hadamard sensing matrices 

Once r and p have been defined for the individual signal and idler subspaces, they may be used 
to construct the corresponding joint space row-selection and permutation vectors, r^j and p^/. 
Consider the construction of rsi first. By Eq. 0, the complete joint space sensing matrix A is 
simply formed by the row-wise Kronecker product of the subspace sensing matrices P 5 and P/. 
As rj and r/ determine the ordering of the Hadamard rows within these patterns, rsi must also 
be a subset of a Kronecker product of and r/. Knowing that the Kronecker product of P 5 and 
P/ will form “blocks” of size [M xN], it is straightforward to show that 

rsi [/] = N (rs [i] - 1) -f r/ \i] (6) 

for i G I...M where r[i] represents the component of r. Note that element-wise counting in 
this article starts at 1 and not 0. 

Because and r/ are chosen at random and will often be over-complete, M > N, and rsi 
will probably have repeating units, and a row within A will appear more than once. This is 
equivalent to taking the same projection more than once, offering no additional information. To 
prevent this, we simply compare each of the values within rsi and eliminate repeating i'* values. 
If rsi[i] is a repeated value, we eliminate r 5 /[/] along with the components r 5 [!] and r/[/]. In this 
way, the number of samples M will decrease yet contain the same amount of information. 

The formation of p^i follows a similar form as r^i, yet it will be of length N^. Although it is 
not a simple Kronecker product, it does follow from the structure in Eq. 0- The structure of 
Ps/ takes the form 

Ps/[A(i-l)-f 7 ] =A(ps[i]-l)-|-p/[;] (7) 

for i G 1...N and j G I...N. Generating randomized Hadamard matrices using r and p for each 
signal, idler, and joint space are summarized below: 

Pj = HAr[rs,P5] 

P/ = H;v[r/,P/] ( 8 ) 

A = H^2[r5/,p5/] 

where the y and x components of H[y,x] refer to the rows and columns of H respectively. The 
construction of A presented in Eq. (j^ allows us to use fast transforms as explained in the next 
section. 

4.3. Joint space fast Hadamard transform operations 

Keeping track of the randomization operations allows the use of fast Hadamard transforms 
when computing A • x and A^ • y. This is accomplished by reordering either x or y according 
to p, taking the fast Hadamard transform, and then picking specific elements from the final 
result according to r. The manner in which they are rearranged and picked depends upon the 
operation A • x or A^ • y in either the data acquisition or reconstruction processes. 

Starting with the data-taking procedure y = A • x, projections are taken in each signal and idler 
system by first constructing individual patterns. Pattern construction is done by fast Hadamard 
transforming basis vectors and then permuting them. Because of the symmetric nature of a 


Hadamard matrix (H = H^), a fast Hadamard transform of a basis vector «[!'], in which the i'* 
component is equal to one and the rest zeros, is equal to the i'* row of the Hadamard matrix 
In short, IK[a[i]] = H[/]. Hence, every pattern P[i,:] can be built according to a fast 
transform by 


H[r[i],:] = 5f[a[r[i]]] 

P[i] = H[r[i],p] (9) 

for i£ 1...M where o:[r[!]], a basis vector whose r[iy^ component is equal to 1, is fast-Hadamard 
transformed and then permuted according to p. 

To take experimental data, many SLM’s, such as digital micromirror devices, are operated 
in a binary fashion of on or off - transmitting light either to or away from a detector. If only 
using one detector per subspace, at any given moment a pattern may only be composed of O’s 
or I’s while Hadamard matrices are composed of I’s and -I’s. To display the full Hadamard 
pattern with one detector per subspace, the data-taking operations must be split into positive and 
negative operations. H ^2 may be decomposed into a sum of four Kronecker products of both 
positive Hat, represented by which is composed of O’s and 1 ’s, and negative H^/, represented 
by which is composed of -I’s and O’s. 

H^2 = (H+ ® H+) + (H^ 0 H^) + (H+ 0 H^) + (H^ ® H+). (10) 

This means that every element in y will require four coincidence measurements. Even though 
4M coincidence measurements are required when using one detector per subspace, the drastic 
sampling performance gained through CS methods is such that 4M ^ N^. Alternatively, if two 
detectors are used in each subspace (one detector to collect projections from the positive com¬ 
ponents and one detector to collect projections from the negative components), the detection 
process could be streamlined to measure each of the four correlations in Eq. ( [TOl l. 

In reconstruction, fast-Hadamard transforms may be utilized by CS reconstruction algo¬ 
rithms to perform the operations A • x and • y. The operation A • x first requires that x be 
inverse-permuted, fast Hadamard transformed, and then have the correct M elements extracted 
from the final result. The inverse-permutation is done by defining an inverse permutation vector 
q as 

q[p[/]]=i (11) 

for all i elements in p. Hence, A • x is realized with the following operations 

y' = 5f[x[q5/]] 

y = y'M- (12) 

The operation A^ • y requires that a vector composed of zeros be filled with the elements 
of y according to r^/, fast Hadamard transformed, and then permuted according to q as follows: 

= y 

x[qs/] = (13) 

Again, these operations work because Hadamard matrices are symmetric. However, it should 
be noted that the true inverse operation is = Hjj/N. When taking the fast transform opera¬ 
tion in Eq. 0, we are explicitly taking the forward fast transform and neglecting the normal¬ 
ization term. Because of this structure, the operations A • x and A^ • y can be utilized by most 
reconstruction algorithms to operate more efficiently. 

The methods outlined in this article can also be applied to joint space signals that are not 
sparse in the “pixel-basis” where 'E ^ 1. Sparse forward and inverse transform operations. 


and need to be applied to x in an appropriate order to bring x and y back into the 

pixel-basis before fast-Hadamard transforming. Hence, the operations A • x and • y become 
A-'I'[x]^' and'I'[A^-y]. 

To reiterate, the novelty presented in this section is in how Eqns. (j^ and Q enable the use 
of fast Hadamard transforms for calculating A • ‘P[x]“^ and 'T[A^ • y] as summarized below. 

. y = A-'T[x]-l 

1. If'T = 1, neglect this step. Otherwise, inverse transform x out of the sparse basis 
using the inverse transform to obtain x' = 'T[x]^*. 

2. Inverse permute x' using q^/ such that x" = x'[q 5 /]. 

3. Fast Hadamard transform x" such that y' = J{[x"]. 

4. Extract M elements from y' using r^/ to obtain y = y^[r 5 /]. 

• x = 'E[A^-y] 

1. Construct a null-vector j3 G 

2. Place the components of y into j3 using r^/ to assign the locations for the elements 
of y such that j3 [r^/] = y. 

3. Fast Hadamard transform j3 such that = TC [)3]. 

4. Permute the elements of P' using psi to obtain x' = j3^[p5/]. 

5. If 'P = 1, neglect this step and let x = x'. Otherwise, transform x' into the sparse 
basis to obtain x = ‘P[x']. 


5. Compressive measurement in a 16.8 x 10^-dimensional correlated space 

5.7. Mutual information 

To demonstrate the practicality of the previous results, we compressively measure and quickly 
reconstruct a 16.8 x 10^ dimensional joint space probability distribution. Up to this point, the 
reason why these joint space measurements are useful has not been explained in detail other 
than to inform the reader of the characterization of correlated systems. When attempting to use 
down-converted photons for information transfer, an important question to ask is. How much 
uncertainty about the position of the signal photon is removed upon knowing the position of the 
idler photon? This quantity is effectively answered by the Shannon mutual information between 
the position statistics of the signal and idler photons I{Xs,Xi) pT) . The mutual information 
quantihes the classical channel capacity and is easily found by hrst measuring the joint space 
probability distribution. Given the discrete random variables’ distributions for signal Xs and 
idler A/, and their allowed set of random values, xs and xj respectively, the mutual information 
is defined as; 



where p{xs,xi) is the joint space probability distribution while p{xs) and p{xi) are the marginal 
probability distributions. In terms of our CS formalism, p{xs,xi) = x. Once x has been procured 
from y, the marginal distributions p{xs) and p{xj) are then found summing over appropriate 
values of p{xs,xj). Since we may approximate the joint distribution as a double Gaussian, we 
may also say 2^(^-5A/) is equal to the Schmidt number of the state which is a measure of the 
number of entangled modes |32 In our case, 2^(^'5A/) is equal to the number of distinct 
channel inputs. 





5.2. Theoretical expectations 


Before reporting our experimental results, it is useful to first estimate the theoretical maximum 
amount of possible mutual information as derived from hrst principles based on the crystal and 
the pump-laser specifications. That result should then be compared with the maximum possible 
information we could measure given our SLM resolution. A thorough calculation characterizing 
degenerate SPDC is done in | [T0| in one transverse spatial dimension. Assuming a double Gaus¬ 
sian bi-photon state, the mutual information in the position domain between down-converted 
photons {Xs and Xj) is 


I{Xs,Xj) = log 2 


97tal+L^,Xp \ 
2c7p 'y/97rZ.^Ap y 


(15) 


where L,, Ap, and Op represent the length of the nonlinear crystal, the pump-laser wavelength, 
and the standard deviation of the Gaussian intensity pump-laser prohle, respectively. We use a 
325 nm pump laser and a 1 mm length nonlinear crystal. The maximum 1/e^ pump diameter 
is listed as 1.2 mm resulting in a Op that is four times smaller, i.e. ffp = 3 x 10 ^ m. The ex¬ 
periment uses approximately degenerate down-converted light because of the paraxial nature of 
beam propagation and the use of narrow-band biters. Our pump laser is approximately Gaus¬ 
sian in two dimensions resulting in a mutual information twice as large as reported in Eq. ( [T5| ). 
We obtain a theoretical maximum mutual information between signal and idler photons of 10.9 
bits. However, when moving from an inbnite-dimensional Hilbert space to a bnite dimensional 
space dictated by the resolution of the SLM and its pixel size, the resulting measurable mutual 
information must be less than or equal to to the continuous variable case of 10.9 bits. 


5.3. Reconstruction algorithm: maximizing the mutual information 

Given that a motivation for these reconstructions is to infer the mutual information shared 
between signal and idler photons, it is reasonable to design a reconstruction algorithm that at¬ 
tempts to maximize the mutual information via the suppression of noise through soft or hard 
thresholding; hard thresholding implies setting components less than a threshold to zero while 
soft thresholding implies brst hard thresholding and then decreasing the amplitude of every 
remaining signal component by the threshold’s amplitude. For practical applications, only the 
largest measurable signal components of x should be used to infer the mutual information be¬ 
cause of the presence of background noise. With this in mind, we use an iterative thresholding 
algorithm fM) that computes the mutual information at each iteration. The program exits when 
the mutual information no longer increases with thresholding. 

The algorithm we use is summarized as follows: 


Xo = c 

Xf+I = fj 2 [xr-{fji [A^-(y-A-x,)]}-|-x,-min(x,)] 


(16) 


where c is a vector composed entirely of ones times a non-zero constant and min(x,) is a vector 
composed entirely of ones times the smallest element in x,. Note that at each iteration we take 
a projection of the current result with the term Tji [A^ • (y — A - x,)]. During the brst iteration, 
A • Xo = 0 because xq = c. Note that A • c = 0 only because we chose to neglect the brst row of 
the Hadamard matrices. The brst iteration results in computing A^ • y as in standard iterative 
style algorithms, ff [...] is an operator that performs soft thresholding on everything within its 


brackets using a biorthogonal 4.4 wavelet transform with a two-level decomposition 135-371. 
Soft threshold in the wavelet basis, often referred to as wavelet shrinkage, is performed using 
the universal threshold of Donoho and Johnstone p8||39l . The bltered signal is then inverse 
transformed back to the pixel basis. 7)2 [...] then performs a hard thresholding on everything 







within its brackets operating in the pixel basis and renormalizes the final result to be a probabil¬ 
ity distribution. The threshold of 772 gradually increases with each iteration. As x, becomes less 
and less noisy through hltering and converges to the true solution, y — A • x, approaches zero. 
However, it never truly reaches zero and results in a small noise term. To prevent injecting ran¬ 
dom noise into each hltered iteration, we take the projection of the noisy term with the current 
clean solution. We then add this projection back to the current solution to prevent discarding 
current signal components. After hard thresholding the hrst iteration, min (x,>o) = 0 where 0 is 
a null vector. 

5.4. Experimental results 

As an experimental demonstration, we compare how the measured mutual information from 
compressive measurements compares to the theoretical mutual information of 10.9 bits. Know¬ 
ing that information from CS resides in the standard deviation of the signal from the different 
projections, this standard deviation must be greater than the shot noise. Otherwise, the signal 
is obscured by the noise. With binary patterns on each SLM we measured coincidence counts 
at a rate of 4 x 10^ counts per second. Since each SLM reduces the incoming flux by ap¬ 
proximately 50%, the total number of coincidences O was approximately four times larger, or 
<!>= 1.6 X 10^ coincidences per second. We chose the integration time such that all four projec¬ 
tions per y-element required a total of 8 seconds due to power constraints. The resulting ratio 
of the standard deviation from the projections relative to the shot noise was 2.4. In addition, the 
sparsity k should be of order N = 4096 due to tight pixel correlations resulting from position 
and momentum correlations, meaning M = 0(Mog(A)) = O(IO^) measurements. We chose 
the number of measurements to be 2 x 10"*^ {M/N'^ « .001) as a reasonable compromise be¬ 
tween the total integration time and reconstruction quality. The resulting scan took just over 44 
hours. To compare these values to a raster scan, the signal-to-noise ratio {SNR) goes as ^J^t/N, 
assuming perfect pixel correlations and uniform illumination. Here, t is the integration time per 
pixel. When raster scanning in an dimensional joint space, the total integration time goes as 
SNR^/<i> for shot noise limited signals. Therefore, a raster scan operating under per¬ 
fect conditions, again considering perfect pixel correlations and uniform illumination, would 
require 50 days to achieve a SNR = 1 for A = 4096. Hence, 50 days represents the bare mini¬ 
mum integration time for raster scans. The reconstructed joint space probability distribution is 
presented in in Fig. 

From the same set of data, a comparison of the recovered mutual information versus the num¬ 
ber of samples M was also conducted. The results state that the resulting mutual information is 
highly dependent on the noise. As the algorithm seeks to hnd the optimal threshold to discern 
the largest number of distinguishable modes, random noise is thresholded into sparse speckle 
patterns when the signal is not large enough to distinguish from the noise. These speckle pat¬ 
terns incorrectly state that there exist tight pixel correlations. To stress the severity of this flaw, 
we consistently recover about 8 bits of mutual information with M —\Q projections. Even for 
large SNR and large M, these speckle patterns may still be present along with the true signal for 
many reconstruction algorithms. Hence, noise will artificially increase the mutual information. 

Reconstruction errors that artihcially increase the mutual information bring into question the 
validity of CS methods for these types of characterizations. However, there exists information 
in the signal and idler marginal probability distributions that can reduce this error. On the right 
hand side of Fig. the marginal distribution for only the signal is shown since the idler dis¬ 
tribution is similar in appearance. The signal’s measured 64 x 64 pixel marginal distribution 
was recovered from the 2 x 10^ projections using photon counts from only that detector. The 
beam is Gaussian; Fig. |^plots everything within the 1 /e^ intensity beam diameter while thresh¬ 
olding the background intensity to zero. With the idler’s marginal distribution taking a similar 
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Fig. 2. Image (a) depicts a zoomed in view that captures the largest components of the 
reconstructed 16.8 x 10® joint space distribution. The largest components are contained 
within the resulting 4-megapixel image shown. The 1 je^ intensity profile, as seen by the 
signal’s SLM, is presented as the signal’s measured marginal distribution in (b). The re¬ 
constructed marginal distribution was obtained through the reconstructed joint space dis¬ 
tribution and is displayed in (c) for comparison. This data was obtained compressively with 
M = 20,000 (.00119 X 64“*) samples in approximately 44 hours and was reconstructed in 
under ten minutes. 


form, these marginal distributions indicate that there are regions where correlations should not 
exist. This information can be used to reduce the error from reconstructions by only allowing 
the algorithm to admit correlations where the marginal distributions say that correlations may 
exist. These results are summarized in Fig. 

It is evident from Fig. [^that background noise significantly alters the results, fabricating a 
larger mutual information. When neglecting marginal information, the recovered mutual infor¬ 
mation values appear to asymptotically decrease as M grows larger, supposedly approaching 
an accurate value. However, including information from the marginals decreases these errors, 
even for small M, by systematically reducing background noise. Note that Fig. plots the 
mutual information in bits. This means that for low sampling percentages when including the 
marginal information, just under 4 bits of mutual information, presumably from noise, is cal¬ 
culated. However, these reconstructions are in a regime where M ^ k\og{N/k) and the signal 
is dominated by noise. 

Using the reconstructed joint space distribution, the marginal probability distributions can be 
calculated and compared to the measured marginal distributions. When not including the meas¬ 
ured marginal information, the joint space distribution becomes riddled with a low-level noise 
that drastically warps the resulting reconstructed marginal distributions. However, including 
the marginal information results in a clean joint space distribution with reconstructed marginal 
distributions that are similar in appearance to the measured values. The reconstructed distribu¬ 
tions in Fig. I^were recovered while including the marginal distribution information. Using 
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Fig. 3. The mutual information, obtained through reconstructions of the joint space dis¬ 
tribution, as a function of the number of projections M is depicted above. The figure lists 
the theoretical maximum mutual information and the corresponding measurement results 
either utilizing or neglecting available information from the marginal distributions. Ne¬ 
glecting information found in the marginals clearly increases the inaccuracy, especially for 
significantly small M, while including this information leads to more realistic values. 


CS and additional information in the marginal distributions, we measure 7.21 bits of mutual 
information, corresponding to 148 correlated modes. The reason we fall short of 10.9 bits of 
measured mutual information is likely due to the low resolution imposed by the beam covering 
a small part of the sample space. This resulted in a discretization of the marginal distributions 
that was not detailed enough to measure all of the available correlations. Another issue is the 
misalignment in the overlap of signal and idler pixels. This causes the appearance of a second 
correlation band as seen in Fig|^(a). Finally, a misalignment in the focus will also result in a 
single pixel in one arm sharing correlations with several pixels in the other arm. 

While we compressively measured the correlations existing between 4096 signal and 4096 
idler SLM pixels, illumination profiles on each SLM suggest that there exist pixels where no 
correlations are even possible. To determine the actual joint space distribution from which we 
could measure possible correlations, it is reasonable to consider the 1 je^ diameter of each Gaus¬ 
sian beam profile and consider the total number of correlations that could exist between those 
pixels alone. The product of these areas, relative to pixel size, then defines the effective joint 
space dimensionality. From these enclosed areas, we calculate a joint space probability dis¬ 
tribution of 3.23 X 10® dimensions. Notice the measured and recovered marginal distributions 
presented in Fig.[^ The marginal distribution recovered from the reconstructed joint space dis¬ 
tribution is smaller than the 1/e^ thresholded distribution displayed immediately above. These 
regions should be identical in the limit of infinite SNR and perfect image-plane optical align¬ 
ment. Again, the difference suggests that either the SNR may not have been sufficient to extract 
all of the existing correlations or that the alignment was imprecise. 


Reconstructions in |21) were limited to a maximum resolution of 10® pixels on a desktop 
computer with 32 GBytes of memory and required several hours to reconstruct x before the 
image quality was sufficient to exit the reconstruction program. However using a laptop with 
8 GBytes of memory, we perform reconstructions of a 16.8 x 10® dimensional joint space in 
under ten minutes with satisfactory results using fast Hadamard transforms with the algorithm 
presented in Eq. (16 1 . We ran the reconstruction algorithm 10 times per sample point to ensure 
























that the results are consistent. The resulting error bars were four orders of magnitude smaller 
than the scale allows and are not shown. It should be noted that the reported mutual information 
values are the result of extracting the largest correlations from the data and are not the result 
of inferring the mutual information from a best Gaussian fit, assuming the resulting SPDC 
follows a commonly approximated double-Gaussian wave function |32 While a Gaussian 
fit would probably characterize the system better and result in a higher mutual information by 
reporting correlations in the tails of the Gaussian, it is not indicative of the actual available 
channel capacity. 


6. Conclusion 

This article describes in detail the methods necessary to efficiently perform high-dimensional 
compressive Kronecker imaging of a joint system. By randomizing Hadamard matrices and uti¬ 
lizing fast Hadamard transforms, CS is performed in the 16.8 million-dimensional Kronecker 
space while containing a 3.23 million-dimensional bi-photon probability distribution. The ex¬ 
periment that would require over 50 days to do a raster scan under perfect conditions is instead 
performed in just under two days and only requires several minutes to reconstruct the data. As 
an example, we compressively measured 7.21 out of 10.9 bits of mutual information while also 
demonstrating how to improve the accuracy of these results utilizing information contained 
within the marginal distributions. We believe these methods will prove to be an invaluable tool 
in measuring the distribution functions of correlated systems as well as other correlation-based 
CS implementations. 
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